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Abstract 



We determine that the sizes of bursts in mean-square current den- 
sity in a reduced magnetohydrodynamic (RMHD) simulation fohow 
a power-law probability density function (PDF). The PDFs for burst 
durations and waiting time between bursts are clearly not exponential 
and could also be power-law. This suffices to distinguish their be- 
haviour from the original Bak et al. sandpile model which had expo- 
nential waiting time PDFs. However, it is not sufficient to distinguish 
between turbulence, some other SOC-like models, and other red noise 
sources. 



1 Introduction 



The widespread occurrence of both self-affine time series with "1//" 
power spectra and spatial fractals in nature led Bak et al. (1987, 1988| ) 
(BTW) to propose the hypothesis of self-organised criticality (SOC) 
( |Bak 19971 ; [Jensen 19981 ; pornette 200C| ) as their common origin. Their 
proposal was based on the demonstration of a "sandpile" cellular au- 



tomaton (see also the earlier work of Katz (1986 )) which appeared to 
be attracted from arbitrary initial conditions ("self-organisation") to 
a critical state characterised by fluctuations on all scales in the energy 
released by the system ("criticality"). Power-law probability density 
functions (PDFs) for the sizes and durations of energy bursts were the 
main observed signatures of criticality, and were tested by finite-size 
scaling of the PDFs with system size ( Cardy 1996| ). 

One of the original applications proposed by BTW for their idea 
was fully developed turbulence, in view of the scaling behaviour of 
such systems, and the intermittency of their energy dissipation. Fur- 
thermore, intermittent turbulence has been an inspiration for later 



"sandpile" -like cellular automata such as the forest fire model of Bak 
et al. (199C). The SOC paradigm has since found many applications 
(|Bak 19971 ; jjensen 1998t [Sornette 2000|) , one of which has been its use 
by Lu and Hamilton (1991| ) and subsequent authors to explain the ob- 



served power-law distributions for the magnitudes, intensities, and du- 
rations of solar flares. SOC has since also been applied to other natural 



and artificial plasma confinement systems, notably the Earth's mag- 
netosphere (a recent review is that of Chapman and Watkins (20011 ) ) 
and tokamaks (e.g. ( [Chapman et al. 2001| )). 

In turn, recent studies in the solar flare context, of both shell 
models ( BofFetta et al. 1999 ) and simulations based more directly on 
the magnetohydrodynamic (MHD) equations ( Georgoulis et al. 1998| ; 
Einaudi and Velli 1999| ) have focused attention on the fundamental 



questions of the ways in which SOC and turbulence may differ (Bof- 



fetta et al. 1999| ) or, conversely, the extent to which SOC may serve as 



a model for turbulence ( Einaudi and Velli 1999| ). These more general 
questions are our main focus in this paper. Also included is some brief 
discussion of the application of these ideas to solar fiares and plasma 
turbulence in the solar wind, magnetosphere, and elsewhere. 

Recent work has shown that magnetically forced 2D MHD turbu- 
lence produces power-law PDFs in the size and duration of bursts in 
spatially averaged Ohmic energy dissipation (r/J^) ( Georgoulis et al. 



"199^ ; [Einaudi and Velh 1999|) (G98). Recall that such power-law PDFs 
are a necessary but not sufficient condition for SOC ( Jensen 199S ). In 
addition Einaudi and Velli (1999[ ) showed that a cellular automaton 
with rules chosen to be consistent with the MHD model also produces 
avalanches (and power-law PDFs). 

In contrast, Poffetta et al. (1999[ ) and [Giufiani et al. (200q ) (B99) 
attributed the presence of a power-law PDF in the observed waiting 
time between solar flares to turbulence. Rather than simply being 
due to the scale-free turbulent cascade itself, B99's suggested mecha- 
nism for the production of power law waiting time PDFs was on-off 
intermittency. They noted that simple prototype models of on-off in- 
termittency and a many-oscillator shell model of turbulence both had 
waiting time PDFs which were power-law, while the original BTW 
sandpile algorithm did not, having instead an exponential PDF of 
waiting times. 

Two important criticisms of B99's model and its interpretation 



have recently been made in the context of solar flares. Einaudi and 



Velli (1999 ) pointed out that its attractors were states in which veloc- 



ity and magnetic field were aligned, whereas they asserted that force 
free (aligned current and magnetic field) states were probably more 
appropriate to the corona. Wheatland (2000) showed how a Poisson 
process varying on a long time scale (e.g. the quasi-periodic solar 
cycle) could convert an exponential waiting time distribution such as 
that from the original sandpile model into a power-law one of the type 
observed. 

These results are however less relevant to the more general plasma 
turbulence case. In physical systems where there is no long-term peri- 



odicity to motivate the "periodic-Poisson" assumption of Wheatland 



"(2000| ) there is no apparent reason to prefer such a mechanism for 
observed power laws to an intrinsically scale-free one. In addition, 
the work of Wheatland (2000 ) seems to have been partly been a re- 
sponse to B99's over-general assertion that all SOC models must have 
exponential waiting time PDFs. This is not true in general (c.f. the 

(2000b[) ; and the 



discussions in Galtier 



(2001[) 



models studied by [Paczuski et al 
interesting questions remain open. 



and Freeman et al. 
1996[) ) 



In consequence, several 



One is whether the burst size and duration PDFs found in simu 
lations of 2D MHD turbulence eoreouhs et al. 19981; lEinaudi andl 



Vehi 19991 ) are also observed for either full MHD or reduced MHD 



(RMHD), in which "slow" kz dependence is retained, (see below). 

The second is whether the power-law waiting time PDFs seen in 
both the B99 one-dimensional shell model and the Galtier (2001) ID 
MHD simulation, are also seen in higher dimensions. 

Finally there is the question of whether a minimal set of scale-free 
"burst" PDFs and power spectra can be identified which suffices to 
identify SOC (or turbulence) in a physical system. This last question 
is also highly topical in magnetospheric and laboratory plasmas (cf. 
Krommes 2000| ; [Freeman et al. 2000b| ; [Kovacs et al. 20011 ). 

In this paper we address these questions by examining time se- 
ries of various spatial averages of the squared electric current density, 
j'^{x,y, z,t), drawn from a RMHD simulation. A threshold method, 
used by previous authors to construct burst-size PDFs, is applied to 
the time series. Power-laws in size (and arguably also in duration) are 
found, extending the forced 2D MHD results of peorgoulis et al. (1998| ) 
and [Einaudi and Velh (1999| ) to forced RMHD. The PDF of wait- 
ing times is also not exponential, confirming that higher-dimensional 
RMHD is in keeping with the predictions of B99 based on a ID shell 
model. Because the fixed threshold method employed detects fractal- 
ity, which was a predicted feature of SOC but is also generic to red 
noise, we then consider to what extent the current evidence is unam- 
biguous. We conclude by suggesting a direction for future research. 



2 Simulation Data and Analysis 

The data analysed here is extracted from a (spectral method) reduced 
MHD simulation which was used in connection with a model for coro- 
nal heating via the coupling of low-frequency Alfven waves and quasi- 
2D turbulence (see [Oughton et al. (2001 ) for further details). Using 



standard (nonlinear) RMHD as a base ( [Montgomery 1982[ ; [Strauss! 



1976; Zank and Matthaeus 199"^ ), the equations were augmented with 



terms representing (i) forcing by a single large-scale Alfvenic mode, 
(ii) reflection of all propagating modes, and (iii) transmission of out- 
ward propagating modes. Physically, one may think of reduced MHD 
as parallel planes of (incompressible) 2D MHD coupled together by a 
strong mean magnetic field (Bq) perpendicular to these planes. Long 
wavelength Alfven waves propagate along the mean field. Thus, the 
fluctuating velocity and magnetic fields (respectively v and b) are 
functions of all three spatial coordinates, but gradients in the Bq di- 
rection are restricted to be weak. Moreover, v and b are strictly 
perpendicular to Bq. 

Here we are primarily interested in various time series characteriz- 
ing such systems. Specifically, those for the spatially averaged mean- 
square electric current density J^(i) = (j^(x, y, z, t))/2, and the k^- 
dependent x- and y-averaged j^/2, denoted as J'^{t,kz), where k^ is 
the Fourier wavenumber in the direction parallel to the mean field 
and angle brackets denote the spatial averaging. Clearly, J^(t) = 



The particular simulation employed has large-scale Reynolds num- 
bers of 800, a resolution of 256^ x 4, (Alfvenic) forcing of the k = 
(1, 1, 1) Fourier mode, and reflection and transmission rates of 1.0 and 
0.3, respectively ( Pughton et al. 2001 ). Note that the reflection and 
transmission parameters are to be interpreted as inverse time scales 
and not as fractions. The simulation was continued for 500 Tb where 
Tb is defined as the time taken for a forced Alfven wave to cross the 
simulation box, which is comparable to the large-scale nonlinear time. 
After a few tens of Tb, the system settles down into a state which is 
more or less statistically steady, reminiscent of the similar behaviour 
of the BTW sandpile (see figure 4.3 of Jensen (1998D ). The time series 
are obtained by calculating the appropriate quantities every 1/10 of 
a box crossing time, after this steady state has been reached. This 
sampling interval was chosen in order to have a manageable amount 
of simulation data. After removing the initial transient, each time 
series J'^{kz,t) consists of 4500 points. The kz = I plane is special, 
since the single forced mode lies in it. 

Figure la shows the statistically steady portion of the time series 
of J'^{t,kz = 2), for which the corresponding PDF and cumulative 
distribution function (CDF) are shown as Figures lb and Ic respec- 
tively. Inspection of the time series itself (Figure la) does not show 
extreme values to the same extent as e.g. figure la of [Georgoulis et aL] 
(1998| ), and so in consequence the PDF we find for it (Figure lb) is 



substantially more symmetric. This is illustrated by the dashed line in 
Figure lb which shows a Gaussian with the same mean and standard 
deviation as those of the time series. 

To measure the distribution of bursts, a fixed threshold method 
was employed, as used by [Freeman et al. (2000b| ) (and previous work- 
ers, see references therein). The size e of a burst was defined as the 
integrated area under the curve between a given upward crossing of a 
fixed threshold and the immediately subsequent downward crossing. 
The duration T was then the time between upward and downward 
crossings, while the waiting time (or inter-burst interval, r) was that 
between a given downward crossing and the next upward crossing. The 
resulting PDF for burst size is plotted in Figure 2a, where the solid 
line indicates the curve corresponding to use of the median value of the 
time series as the threshold, while the eight dotted lines correspond to 
those resulting from the 10th, 20th, . . . , 40th, 60th through 90th per- 
centiles. PDFs for burst duration and waiting time constructed by the 
same method are shown in Figures 2b and 2c. The number of bursts 
bursts thus defined will vary weakly with the threshold chosen, in the 
case of the median threshold the size distribution D{e) was formed 
from 198 events. Despite the symmetrical PDF of shown in Figure 
lb, a power-law PDF is obtained for burst sizes in the range 10~^ 
to 2 units (Figure 2a), which remains stable even as the thresholds 
are varied. Outside this region the points deviate from a power law 
but their statistical weight is low. The dashed line shown is a power 
law fit to those points where the number of samples per bin is greater 
than 4. The PDFs of burst durations and waiting times also resemble 
power-laws in the range 0.2 to 2, beyond which the number of samples 
per bin again falls below 5. Similar plots were obtained for all four k^ 
planes in the simulation, and on averaging over kz- Figure 2d shows 



the PDF for the waiting times plotted on semilog axes, on which an 
exponential distribution would appear as a straight line, confirming 
that the waiting times (as defined herein) for this simulation are not 
exponentially distributed. 



3 Discussion and Conclusions 

We find that the PDF of burst size measured in our simulation is 
power-law in form, independent of the choice of threshold. The PDFs 
of duration and waiting time appear to have the same basic form 
(although the evidence is much less clear cut). Intriguingly, we have 
observed this property in a time series (Figure la) whose PDF (Figure 
lb) is not long-tailed but rather symmetric. Power-law burst PDFs 
are often seen coupled with long-tailed underlying PDFs, and so it is 
sometimes thought that the scaling range of a signal will be controlled 
by the mean to peak ratio (governed by a and /i for Gaussian data, ^ 
for Poisson data, etc.) In fact, however, burst size is also governed by 
the degree of persistence /3 in the signal, because it is not just affected 
by the probability of N points being above a line but is controlled 
by the probability of N non-independent, successive points all being 



above the line. Malamud and Turcotte (1999) have noted that (3 can 



be varied independently of the PDF of the members of a time series, 
so that even a time series with a Gaussian distribution of amplitudes 
but a non-zero /3 would have a non-negligible probability of several 
successive values exceeding a threshold. In addition, for such a series, 
the scaling properties of the distributions of waiting time and inter- 
burst interval are set entirely by {3 and would be power laws. 



This raises an interesting question, however (see also Freeman et al. 
(2000bD ). Red noise ("1//") time series of the type which SOC sys- 
tems were originally expected to produce, and whose appearance SOC 
was proposed to account for, are fractal. The distribution of isosetsQ 
of such a time series is a power-law (Addison 1997). Hence the burst 



duration and waiting time PDFs of red noise when found by a thresh- 
old method should be power-laws, regardless of whether the noise is 
produced by turbulence or one of the class of SOC-type models which 
do produce red noise. The original BTW model can be eliminated as 
its waiting time series was later shown to be uncorrelated and thus not 
red noise ( pensen 1998| ; [Boffetta et al. 1999| ; [Freeman et al. 2000b| ). To 



determine if a process is SOC in the sense of BTW's original proposal, 
one needs information about spatial correlations as well as temporal 
correlations. This is because SOC was proposed as a mechanism link- 
ing spatial fractality and temporal persistence ("l/f noise). This 
reinforces the point that, rather than simply temporal information 
(burst durations and waiting times) or avalanche distributions ( "burst 
sizes"), to differentiate between turbulence and SOC it will be nec- 
essary to make unambiguous predictions about spatial structure for 
each phenomenon, requiring at minimum then availability of spatial 
correlation functions. We note that correlation between bursts e.g. 
the debated "sympathetic" quality of solar flares (Wheatland et al. 



1998| ) could in principle occur through time or space correlation. This 



Defined as the set of times at which the time series crosses a fixed level. 



line of investigation will be pursued in future papers. 
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Figure 1: a) Time series of current density for the reduced MHD simulation. 

b) PDF of current density for the time scries in a). Overplotted is a Gaussian 
distribution with the mean and standard deviation of the time series in Figure 
la. c) CDF of current density for the time series in a). 
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Figure 2: a) Burst size PDF obtained by the threshold method for the time 
series of Figure la). The sohd hue corresponds to the use of the median of 
the time series as the threshold while the dotted lines show the 10th, 20th, 
. . . 90th percentiles. The dashed line passes through those points for which a 
statistically significant number of points is available, b) Burst duration PDF 
obtained by the method of Figure 2a). c) PDF of waiting times between 
bursts obtained by the method of Figure 2a). d) The PDF of figure 2c) 
plotted on a log-linear scale, illustrating that the PDF cannot be fitted well 
by an exponential. 



